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In this paper we present various convective states of zero-Prandtl number Rayleigh- 
Benard convection using direct numerical simulations (DNS) and a 27-mode low- 
dimensional model containing the energetic modes of DNS. The origin of these con- 
vective states have been explained using bifurcation analysis. The system is chaotic at 
the onset itself with three coexisting chaotic attractors that are born at two codimcnsion- 
2 bifurcation points. One of the bifurcation points with a single zero eigenvalue and a 
complex pair (0, ±iw) generates chaotic attractors and associated periodic, quasiperi- 
odic, and phase-locked states that are related to the wavy rolls observed in experiments 
and simulations. The frequency of the wavy rolls are in general agreement with u of 
the above eigenvalue of the stability matrix. The other bifurcation point with a double 
zero eigenvalue produces the other set of chaotic attractors and ordered states such as 
squares, asymmetric squares, oscillating asymmetric squares, relaxation oscilla tions with 
interm ediate squares, some of which are common to the 13-mode model of 



Pal et al 



( 20091) 



Key Words: Buoyancy-driven instability, Low-dimensional models, Bifurca- 
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1. Introduction 



Origin of instabilities, patterns, and chaos in Rayleigh-Benard convection (RBC) is an 



interesting problem in fluid dynamics ( Chandrasckhar 



1987 



Manneville 



2004; 



Ahlers. Grossman fc Lohse 



1961 



Busse 



1981 



Bhattachariee 



2009). Prandtl number, P (ratio of 



kinematic viscosity v and thermal diffusivity k) and Rayleigh number, R (ratio of buoy- 
ancy and dissipative terms) are the two critical parameters for RBC. Some of the key 
and difficult problem in this field are related to the onset of convection for zero-Prandtl 
number (zero-P) and low-Prandtl number (low-P) fluids. For zero-P and low-P convec- 
tion, the nonlinear term of the Navier-Stokcs equation plays an important role, and it 



secondary oscillatory patterns (jBusse 



Busse 


1972 


Clever & Busse 


1974 


Busse & Bolton 


1984) 



In the present paper we will study various patterns and chaos for zero-Prandtl num- 
ber convection using bifurcation analysis. We use direct numerical simulations and a 
low-dimensional model for this purpose. 

Low-Prandtl number fluids, for example, mercury (P w 0.02), liquid sodium (P rj 



terns and chaos ( 


Croquette 


1989a 1 


] ; Cross & Hohenbers 


1993 


Bodenschatz. Pesch & Ahlers 


2000 


Das. Ghosal & Kumai 


2000; 


Stein & NorlunJ 


2000 


). In terrestrial experiments, the 



Prandtl number cannot go below that of liquid sodium (~ 0.01). In addition, the visu- 
alization of flow patterns inside low-P fluids like mercury and sodium is quite difficult. 
Due to these limitations of experiments, numerical simulation s are v e ry sig nificant in the 



study of low-P and zero-P convection. The numerical work of 



Thuall (jl992f) indicate that 



the properties of low-P convection as P — > are quite close to that of zero-P convection. 
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Hence low-P convection as P — > appears to approach zero-P convection as a limiting 

case. Therefore, zero-P convection is very useful for understanding the properties of low- 
P convection. Even though numerical analysis of zero-P convection is quite tricky due 
to its inherent instabilities near the onset, it provides certain advantages. For P = the 
thermal modes are slaved to the velocity modes, hence the number of independent vari- 
ables are less than those required for low-P analysis. Also, the time steps for numerical 
simulations of very low-P convection could be very s mall due to th e stiffness of the equa- 



1992). These features led us 



tions, which is not a limitation for zero-P convection (|Thuall l~ 
to investigate zero-P convection in detail for understanding the convective patterns and 
chaos in low-P fluids. 



Thuall (|l992l ) was one of the first to simulate zero-P convection for free-slip and no-slip 
boundary conditions. He reported various supercritical oscillatory instabilities, regular 
and chaotic patterns, etc. in his simulations. The patterns observed by Thual are two- 
dimensional rolls, periodic and qua siperiodic rolls, squares, travelling waves, etc. For 



low-P fluids, 



Meneguzzi et al. 



(1987|) performed numerical simulations for P — 0.2 under 
stress-free boundary conditio ns, and for P = 0.025 under no-slip boundary conditions, 
and observed similar patterns. IClever &: Bussd (|1990i ) found travelling wave solutions for 



low-P convection under no-slip boundary conditions. These observations indicate that 
zero-P convection contains pertinent features of low-P convection. Also, free-slip and 
no-slip boundary conditions exhibit similar convective flow patterns. 

The oscillatory instabilities, and regular and chaotic patterns found in numerical 



and liquid sodium ( 


Rossbv 


969; 


•Crishnamurti 


1970 




1973 




Maurer & Libchaber 


1980 


Libchaber. Laroche & Fauvt 


1982 


Libchaber. Fauve & Laroche 


1983 


) . Some of the most 



commonly observed patterns in experiments are stationary, periodic, quasiperiodic, chaotic, 
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and travelling rolls, as well as squares, asymmetric squares, and phase locked states. 
Chaotic rolls have been observed to appear through both period doubling and quasiperi- 
odic routes in some of these experiments- 



Con vection experiments of 



Willis fe Dcardorff 



199(1 ) 



1974 ) , and numerical simulations of 



Meneguzzi et al 



(1987), 



Lipps 



Thual 



1976), 



19921 ). 



(19G7) 



Krishnamurti 



Bolton fc Busse 



Ozoe&Hara 



1970 



19851), 



1973), 



Busse fe Whitehead 



Clever fc Busse 



1995) indicate the presence 



of oscillatory instability of the two- dimensional convective rolls and the resultant sta- 



ble w avy rolls. Busse and coworkers (jBusse 



1972 



Clever fc Busse 



1974: 



Busse fc Bolton 



19841) showed using perturbative analysis that the two-dimensional rolls become unsta- 



ble to oscillatory three-dimensional disturbances w hen the ampl itude of the convective 



motion exceeds a finite critical value. According to 
instability takes the form 



Bussej ( 19721 ). the condition for the 



Rt 

Rr 



1 « 0.310P 2 



(1.1) 



where R c is the critical Rayleigh number for the onset of convection, and Rt is the 
Rayleigh number where oscillatory instability starts. In addition, the time period of 
the oscillati ons measure d in the units of d 2 /v (viscous time scale) is indep endent of P 



(Eq. (5.2) of 



Busse 



19721 ) . In a related work. 



Fauve. Bolton fc Bracheti (jl987l ) investigated 



the origin of instabilities in low-P convection using phase dynamical equations and argued 
that the instability always saturates into travelling waves. In the present paper we will 
explore the origin of oscillatory instability using numerical simulations and bifurcation 
theory. 

Origin of plethora of convective patterns observed in convective experiments and nu- 
merical simulations can be quite intricate. Each run of a convective simulation takes 
significantly long time, so it is not possible to scan the parameter space minutely for 
deciphering detailed bifurcation scenario. Experiments too have their complexities and 
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limitations. The large number of modes present in simulations and experiments tend 

to obscure the underlying dynamics. These difficulties are circumvented by a powerful 
and complimentary approach in which the system is analyzed using an appropriately 
constructed low-dimensional model. Relatively low computational cost for running low- 
dimensional models, and the ease of construction of the bifurcation diagrams are some 
of the distinct advantages of these models compared to the experiments and simulations. 
However care must be taken to take all the relevant modes of the system for constructing 
the low-dimensional models. 

Number of active modes near the onset of convection is not very large, so low- dimensional 
model s consisting of these active modes are very useful for analyzing this regime. 



Kumar. Fauve fc Thual 



( 1996) showed using a 6- mode model of zero-P convection that the growth of the 2D rolls 



satura te through the generation of vertical vorticity (wavy nature). 



Kumar. Pal fc Fauve 



( 2006) observed critical bursting in the above model during the saturation. 



Pal fc Kumar 



20021 ) explained the mechanism of selection of the square patterns using a 15-mode 



model of the zero-P RBC. Recently 



Pal et al. 



(|2009h constructed a 13-mode model us- 



ing the energetic modes of DNS, and performed a bifur cation ana l ysis u sing the model 



Pal et al. 



20091 ) could explain 



and simulation results. Using the bifurcation diagram, 
the origin of squares (SQ), asymmetric squares (ASQ), oscillating asymmetric squares 
(OASQ), relaxation oscillations with an intermediate square regime (SQOR), and three 
kinds of chaotic attr actors. In ad dition, the above patterns were observed in both DNS 
and models. Earlier iThual J 19921) had observed SQ and SQOR patterns in his simula- 

1984 ) studied the transition from 2 D rolls to square patterns 



tions. 



Jenkins fc Proctoi 



Knoblochl (jl992f ) studied 



in RBC for different Prandtl numbers using analytical tools, 
the stability of the SQ patterns using a complementary procedure called the amplitude 
equations. 
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Pal et al. 



K. Vnr ma 

( 20091 ) was the absence of wavy 



A major limitation of the 13-mode model of 
rolls near the onset. This limitation is overcome by extending this model to a 27-mode 
model by incorporating the corresponding dominant modes. In the present paper we 
perform a bifurcation analysis of this model and explore the origin of the various con- 
vective patterns of zero-P convection with special emphasis on the oscillatory instability 
and related wavy roll patterns. All the features of the 13-mode model are reproduced 
in the 27-mode model by construction. We will show in our discussion that properties 
of the wavy rolls of the 27-mode model matches reasonably well with those observed in 
experiments and simulations. 

The organization of the paper is as follows: In section [5J we describe the basic hydro- 
dynamic system considered for the study. The low dimensional model will be derived 
in section [3l Section |4] contains the results of the numerical simulations and the low- 
dimensional model. Various bifurcation diagrams are described in this section. A brief 
study of the wavy rolls observed in the low-dimensional model is presented in section [5l 
We finally conclude in section El 



2. Governing equations and direct numerical simulations 

The RBC system consists of a conducting fluid of kinematic viscosity v, thermal diffu- 
sivity k, and coefficient of volume expansion a confined between two conducting plates 
separated by a distance d and heated from below. In the zero-P limit the equations under 



Boussinesq approximation take the form ([Spiegel 



1962 



Thual 



19921) 



d t (V 2 v 3 ) = V 4 v 3 + RV 2 H 9 - e 3 • V x [(w-V)v - (v-V)w] , (2.1) 
d t uj 3 = V 2 w 3 + [(«-V)«3 ~ (v-V)w 3 ] , (2.2) 
V 2 e = -V3, (2.3) 
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V-v = 0, (2.4) 

where v(x,y, z,t) = (vi,V2,V3) is the velocity field, 9(x,y,z,t) is the deviation in the 
temperature field from the steady conduction profile, u> = (u>i, U2, W3) = V x v is the 
vorticity field, §3 is the vertically directed unit vector, and = d xx + d vy is the horizon- 
tal Laplacian. The equations have been nondimensionalized using d as the length scale, 
d 2 jv as the time scale, and vftd/ n as the temperature scale, where j3 is the uniform tem- 
perature gradient. The two nondimensional parameters in the equations are the Rayleigh 
number R = af3gd 4 / 'vk and the Prandtl number P = v j k, where g is the acceleration 
due to gravity. In the following discussions we will also use the reduced Rayleigh number 
r = R/ R c as a parameter. 

We consider perfectly conducting boundary conditions for the top and bottom plates 
along with the free-slip boundary condition for the velocity field. Consequently 

v 3 = d s vi = d 3 v 2 = 9 = 0, at z = 0, 1. (2.5) 

We assume periodic boundary conditions along the horizontal directions. 

Equations (|2.HI2.4j) are numerically solved using direct numerical simulations (DNS) 
under the above boundar y conditions (Eqs. 12 51). DNS wer e performed using a pseudo- 



spectral code TARANG (jCanuto et al 



1988; 



Verma 



2008h in a box with aspect ratio 



= 2y/2,T y = 2^2. Various grid resolutions, 32 x 32 x 32, 64 x 64 x 64, have been 
used for the simulations. These grids provide well-resolved simulations near the onset 
of convection. We used fourth-order Runge-Kutta scheme (RK4) for time advancement. 
Each run was carried out till the system reaches a steady state. The DNS runs were 
performed for 0.98 < r < 1.25. 

We perform around 200 simulation runs, yet it is insufficient to construct the bifurca- 
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tion diagram of the system in detail. For this purpose we construct a low-dimensional 
model, which will be described in the next section. 

3. Low-dimensional model 

We identify 27 Fourier modes that have significant energy (approximately 1% or more 
of the total energy) in the DNS runs and create a low-dimensional model. The modes of 
the model are shown in Fig.[T]and they account for approximately 98% of the total energy. 
The triangles represent the interacting triads. Note that only some of the interacting 
triads of the model have been shown in the figure. Care has been taken to include sufficient 
number of modes so that the model reproduces many features found in experiments and 
simulations. Also, we ensure that the model results are reasonably close to the simulation 
results. The vertical velocity field (V3) and vertical vorticity field (W3) in terms of the 
chosen modes are 

v 3 = Wioi(t) cos(kx) sin(7T2;) + Wmi(t) cos(ky) sin(7rz) 

+W 20 2(t) cos(2fcx) sin(27rz) + W022CO cos(2fcy) sin(27rz) 

+Wios(t) cos(kx) sin(37rz) + Woi3(t) cos(ky) sin(37rz) 

+W / 3oi(i) cos(3A;a:) sin(7rz) + Wo3i(t) cos(3fcy) sin(7rz) 

+Wi2i(t) cos(kx) cos(2ky) sin(7rz) + W2ii(t) cos(2fca;) cos(ky) sin(7rz) 

+Wn2(t) cos(kx) cos(ky) sin(27rz) + Wm(t) sin(fcx) sin(fcy) sin(7rz) (3.1) 

^3 = Z wo (t) cos(kx) + Z QW cos(fcy) 

+Zno(t) sin(fcx) sin(fcy) + Z\n(t) sin(fcx) sin(fcy) cos(27rz) 
+Z3io(t) sin(3fcx) sin(fcy) + Z\^{t) sin(fca;) sin(3fcy) 
+Z120 (t) cos(fcx) cos(2fcy) + Z2io(t) cos(2fca;) cos(fcy) 
+Z W 2(t) cos(kx) cos(27rz) + Z i2(t) cos(ky) cos(27rz) 
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+^20i(t) cos(2fcx) cos(7rz) + Zo2i(t) cos(2ky) cos(ttz) 

+Zm(t) cos(kx) cos(ky) sin(7rz) + Z 12 i (i) sin(fca;) sin(2fcy) cos(7rz) 

+Z 2 ii(t) sin(2fcx) sin(fcy) cos(7rz) (3-2) 

where Wi mn and Zi mn are the Fourier amplitudes of the vertical velocity and vertical 
vorticity modes respectively with the three subscripts (/, m, n) indicating the wavenumber 
components along the x, y, and z directions respectively. The modes (1, 0, 1) and (0, 1, 1) 
are the most important modes of our model, and they represent the rolls along y and 
x directions respectively. For the square pattern, the most important participating triad 
is {(1, 0, 1), (0, 1, 1), (1, 1, 2)}. Note that the wavenumbers of the interacting triad satisfy 
k = p + q. 

The horizontal components of the velocity field can be computed using the incom- 
pressibility condition of the velocity field (Eq. (j2.4[l ). and the temperature field 9 can 
be computed using Eq. f|2.3|) . A Galerkin projection of the RBC equations (|2. 1112. 21) on 
the above modes provides a set of 27-dimensional coupled first-order nonlinear ordinary 
differential equations for the amplitudes of the above Fourier modes. On this 27-mode 
model, we perform a detailed bifurcation analysis that we describe in the subsequent 
sections. 

Our 27-mode model is a superset of the 13-mode model of 



Pal et al. 



(2009). The 



nonlinear interactions of the 13-mode model are indicated by blue curves in Fig. [TJ 
Additional interactions induced by the new modes of the 27-mode model are represented 
by red curves in the figure. A primary motivation for the 27-mode model is to be able to 
generate wavy rolls. The triads {(0, 1, 0), (1, 0, 1), (1, 1, 1)} and {(1, 0, 0), (0, 1, 1), (1, 1, 1)} 
play a critical role in inducing wavy rolls along the y and x axes respectively. In the 
present paper we will investigate the dynamics of these wavy rolls in RBC using numerical 



10 



P. Pal, S. Paul, P. Wahi and M. K. Verma 

(1,1,0) 



(2,0,1) 



(0,2,1) 



(0,1,3> 




(1,0,3) 



(3,0,1) 



(0,3,1) 



(3,1,0) (1,3,0) 

Figure 1. The interacting modes of the 27- mode model. Interactions corresponding to the 
13-mode model are shown in blue, whereas new interactions unique to 27-mode model are 
shown in red. The modes in green are active in the wavy roll flow patterns. 

simulations and the 27-mode model along with other features that are generated by the 
inclusion of these modes. 

We numerically solve the 27-mode model by employing accurate ODE solvers of MAT- 
LAB. As a result we observe a variety of convective patterns: squares (SQ), asymmetric 
squares (ASQ), oscillating asymmetric squares (OASQ), relaxation oscillations with an 
intermediate square regime (SQOR), wavy rolls, chaotic squares, etc. We have illustrated 
three snapshots each of OASQ in Fig. [21 SQOR in Fig. 02 and wavy rolls in Fig. |U For 
dynamics of these patterns as well as other patterns mentioned in this paper refer to 
the acc ompan y ing vi deos. Note that all the above patterns were also found in our DNS. 



Earlier 



Thuall (jl992l ) in his DNS of zero-P convection had shown the existence of SQ, 
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(a) (b) (°) 




0.2 0.4 0.6 0.8 1 



Figure 2. Oscillating asymmetric square (OASQ) pattern in the mid-plane of the convection 
box in zero-P RBC. The pattern is obtained from the 27-mode model at r = 1.1. Snapshots at: 
(a) t = 0, (b) t = T/4, and (c) t = T/2, where T is the time period of oscillation. 



(a) (b) (c) 




Figure 3. Relaxation oscillations with an intermediate square regime (SQOR) pattern 
observed in the 27-mode model at r — 1.05. Snapshots at: (a) t = 0, (b) t = T/4, and (c) t = 

T/2. 

SQOR, oscillatory quasihexagons (SQOS), chaotic squares (SQCH), and chaotic quasi- 
hexagon (HXCH). Thual observed the oscillatory and chaotic quasi-hexagons for Rayleigh 
numbers beyond the range investigated in this paper. 

We investigate the origin of various convective flow patterns from the bifurcation dia- 
grams generated using the low dimensional model. To generate the bifurcation diagram, 
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0.2 0.4 0.6 0.8 



Figure 4. Wavy roll pattern in the 27-mode model at r — 1.15. 
Snapshots at: (a) t = 0, (b) t = T/4, and (c) t = T/2. 

we start first with the fixed points of the system. We compute the fixed points using 
the Newton-Raphson method for a given r, and these fixed points are subsequently 
continu ed using a fixed arc-length based continuation scheme fo r the neighbouring r 



values (jNandakumar fc Chatteriee 



2005 



Wahi fc Chatteriee 



2008(1 . The stability of the 



fixed points are ascertained through an eigenvalue analysis of the Jacobian. New branches 
of fixed points and limit cycles are born when the eigenvalue(s) become zero (pitchfork) 
and purely imaginary (Hopf) respectively. This process is continued on the new branch. 
For aperiodic and chaotic solutions, we resort to numerical integration and report the 
extremum va lues of the important modes. We use our own MATLAB code as well as 



MATCONT (|Dhooge. Govaerts fc Kuznetsov 



2003) for the analysis. 



4. Bifurcation analysis using model and simulation results 

In the present section we numerically solve Eqs. (|2. 1112.41) using DNS and the 27-mode 
model in the range 0.98 < r < 1.25. This range of r values is near the onset of convection. 
We will present the bifurcation diagrams associated with the different attractors using 
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Figure 5. Three dimensional view of the bifurcation diagram exhibiting the fixed points only. 
Solid and dashed curves represent the stable and unstable fixed points respectively. Black, 
blue, and cyan curves represent stationary squares (SQ), asymmetric stationary squares (ASQ), 
and conduction state respectively. Solid purple lines on the VK101 and M^on axes represent the 
neutrally stable 2D roll solutions and chained purple lines represent unstable 2D rolls. CBf 
(black diamond), CB2 (red dots) and CB3 (green squares) are the codimension-2 bifurcation 
points on the r = 1 plane. The critical eigenvalues at CB1, CB2 and CB3 are denoted by A c in 
the figure. 

the low-dimensional model followed by a detailed comparison of the model results with 
those obtained from DNS. 
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4.1. Fixed points of the system 

Fixed points of a dynamical system and their bifurcations provide important clues about 
the system dynamics. Therefore, we start our analysis by locating all the fixed points 
of the 27-mode model. In Fig. [5] we display the projection of these fixed points on the 
Wioi-Won plane as a function of r. For r < 1, the only stable fixed point of the system 
is the origin, which corresponds to the pure conduction state (cyan curve). At r — 1, the 
conduction state loses stability, and neutrally stable pure roll solutions (purple curves of 
Fig. [5]) and four unstable branches of symmetric square solutions (dotted black curves of 
Fig. [5]) satisfying |Wioi| = |Woii| are born. Note that the stability matrix at the origin 
(CB1) has double zero eigenvalues, and hence CB1 is a codimension-2 bifurcation point. 

The neutrally stable 2D roll solutions of the r = 1 plane become unstable through 
another codimension-2 bifurcation at CB2 {(W W i ~ ±13.44, Won = 0); (Won ~ ±13.44, 
Wioi =0)}, which is represented by large red dots on the r = 1 plane. The stability 
matrix at CB2 has a zero eigenvalue and a purely imaginary pair (A c = (0,±iu>) with 
u) ss 14.2). As a consequence of the complex eigenvalues of CB2, for r > 1 periodic 
solutions are born, and the 2D rolls lose their stability (chained purple lines of Fig. [5]). 
These periodic solutions are also unstable akin to the symmetric square fixed points 
associated with CB1. We will show later that the wavy rolls are associated with CB2. 



Note that CB2 and their associated attractors arc absent in the 13-mode model (Pal et 



2009) 



The unstable roll solutions which persist at r = 1 with amplitudes greater than 
13.44 subsequently undergo yet another codimension-2 bifurcation at CB3 {(Wioi = 
±26.94, Won = 0); (Wioi = 0, Won = ±26.94)}, which are represented as large green 
squares on the r — 1 plane. The stability matrix at CB3 has double zero eigenvalues (A c = 
0,0). Consequently unstable asymmetric square solutions with |Wioi| 7^ |Won| (dotted 
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blue curves) are born. This bifurcation is also present in the 13-mode model (|Pal et 



20091) . and the attractors associated with this bifurcation in the 13-mode model carry 
over to the 27-mode model as well. Note that the ASQ solutions of the 27-mode model 
have two pairs of unstable eigenvalues as opposed to one unstable pair in the 13-mode 

model due to the prese nce of CB2. 

For low-P convection iMishra. Wahi fc Vermal <|2010f ) have earlier observed that the 2D 



rolls undergo a pitchfork bifurcation followed by a Hopf. In the limiting case of zero-P, 
the Hopf bifurcation point merges with the pitchfork bifurcation point at r — 1, and 
the critical eigenvalues at this bifurcation point are (0, 0) giving rise to a codimension-2 
bifurcation CB3. With an increase in r, the double zero eigenvalues split into unstable 
complex co njugate pair giving rise to a scenario very simil ar to the Takens-Bogdanov 



bifurcation ( Guckenheimer fc Holmes 



1983 



Kuznetsov 



1995). 



As shown in Fig. [SJ on the projection of W101-W011 plane, there are 13 unstable 
fixed points for r > 1: one corresponding to the pure conduction state, four satisfy- 
ing |Wioi| = | Won | (SQ), and the remaining eight satisfying |Wioi| ^ |Won| (ASQ). 
After two successive inverse Hopf bifurcations, to be described later, the unstable ASQ 
fixed points become stable. These stable fixed points are shown by the solid blue curves in 
Fig. [5] Subsequently at r ~ 1.1690, these stable ASQs merge (via a pitchfork bifurcation) 
with the symmetric square solutions that originate from CB1 and stabilize them (solid 
black curves in Fig. [5]) . Note that there is a small difference in the values of r correspond- 
ing to the stabilization of the ASQ solutions and SQ solutions for the 13-mode model 
and the 27-mode model. In summary, the 3D figure (Fig .EQ) is qualitatively similar to the 



2009) except the codimension-2 



corresponding figure of the 13 mode model (|Pal et al. 
bifurcation at CB2 which is responsible for the wavy rolls. In the following discussions 
we will describe the bifurcation diagrams including limit cycles, chaotic attractors etc. 
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1.05 



1.15 



1.25 



Figu re 6. Bifurcati on diagram Bif-13M of the 13-mode model for 0.98 < r < 1.25 (adapted 
from Pal et al 



2009 ). The stable branches corresponding to SQ and ASQ are represented by 
solid black and blue lines respectively. Red, green, and brown points represent the extrema of 
OASQ, SQOR, and chaotic solutions respectively. A zoomed view of the bifurcation diagram for 
the chaotic regime is shown in the inset. Branches corresponding to the unstable fixed points 
are represented by dashed lines. Cyan line represents the conduction state. 

Bifurcation diagrams of the 27-mode model are quite complex. They include six differ- 
ent types of chaotic attractors, various types of fixed points and periodic solutions. The 
model also has multiple coexisting attractors for a given value of r. To disentangle its 
complexity, we present the bifurcation diagrams as four separate diagrams, "Bif-13M", 
"Bif-A", "Bif-B", "Bif-C", that highlight different features of the dynamics. First we 



draw the bifurcation diagram "Bif-13M" associated with the 13-mode model (jPal et al 



20091 ) . which is a subset of the 27-mode model. Later, we will contrast the bifurcation 
diagram of the 27-mode model (denoted by "Bif-A") with the diagram of the 13-mode 
model. 
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4.2. Bifurcation diagram of the 13-mode model 
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As illustrated by Fig. [TJ the 13-mode model is a subset of the 27-mode model. If we force 
only the modes of the 13-mode model to be nonzero, and others to be zero, naturally 



the bi furcation diagram corresponding to the 13-mode model is reproduced. 



Pal et al 



(|2009r ) contains a detailed discussion on this diagram and the associated flow patterns 
(both from DNS and the model) . Here we only provide a brief description. 

The bifurcation diagram for the 13-mode model is illustrated in Fig. [6] in which we 
plot the positive value of \Wioi\ extre mum as a function of r. For zero-P convection, chaos 



is observed at the onset itself. Therefore, 



Pal et al. 



([20091 ) start their analysis at r = 1.4 



where symmetric square (SQ) patterns are observed. These states are represented by the 
solid black curve of Fig. [6] (here the diagram is shown only for r < 1.25). At around 
r s» 1.2201, SQ branches bifurcate to ASQ solutions through a supercritical pitchfork bi- 
furcation. The solid blue curves of Fig. |6] correspond to ASQ patterns. The ASQ branches 
bifurcate to OASQ solutions (the solid red curves) through a Hopf bifurcation. The limit 
cycles thus generated grow in size and touch the saddles (dashed line of the SQ branch) 
to create a very narrow window of homoclinic chaos. After this, the system again becomes 
periodic (SQOR) with the merger of the limit cycles. The SQOR patterns transform to 
a chaotic attractor Chi through a homocl i nic bi furcation. Chi turns to Ch2 and subse- 



quently to Ch3 th rough "crisis" . 



Pal et al. 



and DNS. Earlier, 



([20091 ) observed these patterns in both model 



Thuall (|1992l) had observed SQ, ASQ, and SQOR patterns in his DNS 



runs. 



In the subsequent subsections we will describe the bifurcation scenario for the 27-modc 
model. 
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4.3. Bifurcation diagram Bif-A of the 27-mode model 

The square pattern described above is also observed in the 27-mode model. However, for 
SQ in the 27-mode model, twelve modes Win, Z 1W , Z\yz, Z im , Z$ w , Zm, Z 2 io, Z V2 q, 
Z201, Z021, Z102, and Z012 still remain zero. When we continue the SQ branch of the 27- 
mode model, we obtain a new bifurcation-diagram called "Bif-A" shown in Fig. [7] The 
bifurcation diagram Bif-A is qualitatively similar to Bif-13M except in a narrow window 
of 1.116 < r < 1.128 where additional bifurcations are observed. The ASQ branch of 
solutions in Bif-A have seventeen active modes with the modes Wm, Z W q, Zqiq, Zm, 
Z210, Z120, Z 2 qi, Z 2i, Z W 2, and Z i 2 as zeros. 

The new features of Bif-A are as follows. At r = 1.1260, the ASQ branch undergoes 
a supercritical Hopf bifurcation (HI, see Fig. [8]) resulting in a time-periodic convective 
flow as illustrated in Fig. |9^a,b), where we show a projection of the limit cycle obtained 
from the DNS and the model on the Wm-Zoio plane. All the 27 modes are active for 
these periodic flow patterns. 

As r is reduced further, at r = 1.1257 a new frequency incommensurate with the 
original frequency is born through a Neimark-Sacker bifurcation (NS1) and the limit 
cycle becomes unstable. Here, a pair of imaginary Floquet multipliers cross the unit 
circle outwards as illustrated in Fig. irUT a 1 ). The phase space trajectory of the system 
on the VFm-Zoio plane is therefore quasiperiodic as demonstrated in Fig. Etc,d) for the 
DNS and the model respectively. The unstable limit cycle continues till r — 1.0651 where 
they meet the unstable limit cycles of the wavy rolls, which will be discussed in § 14.51 

On further reduction of r, at r = 1.1226, another Hopf bifurcation (H2) takes place 
on the unstable ASQ branch. The limit cycle born from H2 is however unstable. At r = 
1.1181, this limit cycle becomes stable via an inverse Neimark-Sacker bifurcation (NS2) 
wherein the unstable Floquet multiplier pair enters the unit circle as evident in Fig. llOf b). 
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Figure 7. Bifurcation diagram Bif-A of the 27-mode model with the same color convention as 
Bif-13M shown in Fig.[S] This diagram is qualitatively similar to Bif-13M. New features of Bif-A 
are the HI, NS1, H2, and NS2 bifurcations shown in the boxed region (1.116 < r < 1.128) whose 
zoomed view is shown in the inset. 
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Figure 8. Plot of Win vs. r near the first Hopf bifurcation (HI) of the ASQ branch. The 
solid brown curve represents the limit cycles generated after the first Hopf (HI). 
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Figure 9. Projection of the phase space on the W\xi - Zoio plane. The limit cycles born through 
HI in (a) DNS at r=1.138 and (b) model at r=1.1258. Quasiperiodic attractor born through 
NS1 in (c) DNS at r = 1.131 and (d) model at r = 1.1245. 




Re(Floquet Multiplier) 



Figure 10. Movement of Floquet multipliers in the complex plane during NS1 and NS2 bifur- 
cations, (a) At NS1, a pair of complex Floquet multipliers move out of the unit circle (red dots 
changing to blue), (b) At NS2, a pair of complex Floquet multipliers enter the unit circle (red 
dots again becoming blue). 
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Figure 11. Projections of the two coexisting OASQ solutions on the Wioi -Won plane just 
after the NS2 bifurcation point: (a) in DNS at r — 1.0768 and (b) in the model at r = 1.095. 

The resulting stable limit cycle is the oscillatory asymmetric square (OASQ) solution of 
the 13- mode model. Figure QT] shows a projection of the limit cycle corresponding to the 
OASQ solution obtained from the DNS and the model on the Wioi-Woii plane. The 
quasiperiodic solutions exist only in the range r = 1.1181-1.1257, i.e., between NS1 and 
NS2 and disappears after NS2. Note that the attractors between HI and NS2 contain all 
the 27 modes, but beyond NS2 the ten modes Win, -Zioo: Zoio, Zm, Z210, ^120, -Z021, 
Z102, and Zqi2 again become zero. 

Beyond NS2, the patterns and associated bifurcations from OASQ to SQOR to the 
chaotic attractors Chi, Ch2 and ChS in the decreasing r direction are exactly the same as 
for Bif-13M. The range of r corresponding to these patterns are approximately r = 1.086- 
1.1181 for OASQ, r = 1.0046-1.086 for SQOR, r = 1.0034-1.0046 for Chi, r = 1.0025- 
1.0034 for Ch2 and r = 1-1.0025 for Ch3. Comparison with the 13-mode model indicates 
that the ranges of r for the above patterns as well as that for SQ and ASQ arc different. 
This is due to the fact that more than 13 modes are active in the present model for 
Bif-A. The 27-mode model reproduces the ranges of the patterns obtained from DNS 
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more accurately. Note that the modes absent for SQ, ASQ, OASQ, SQOR, Chi to ChS 
flow patterns are active for wavy rolls that will be described in i)4.5l 

4.4. Bifurcation diagram Bif-B of the 27-mode model 

When we start from an arbitrary initial condition for r > 1 near the onset, most often 
the system tends to another attractor Ch4, which differs from Chi, Ch2, and ChS. The 
codimension-2 bifurcation point CB3 of Fig. [5] generates the chaotic attractor Ch4 as 
well. The bifurcation diagram "Bif-B" shown in Fig. [TS] contains the attractor Ch4- The 
attractor Ch4 coexists with the chaotic attractors Chi, Ch2, Ch3 and SQOR for r = 1- 
1.056. This feature is illustrated in Fig. [T31 where a phase space projection (both from 
the model and DNS) for two different initial conditions at r = 1.0342 yield the SQOR 
(green curve) and the Ch4 (grey trajectory) attractors. Clearly the trajectories of Ch4 
explore all the four quadrants of the W101-W011 plane. The qualitative behavior of the 
Ch4 attractor is similar to Ch2, but its size is larger than Ch2 (compare with Fig. 5 of 



Pal et al 



2009) 



Note that there are two complex pairs of unstable eigenvalues associated with ASQ. The 
first one carries from CB2, and the second one is due to the splitting of the double zero 
eigenvalues at CB3 into a complex pair (see §4.ip . The Ch4 attractor is associated with 
the first pair, while the ChS attractor is related to the second pai r. The disappe arance of 



19941 ) wherein 



the chaotic attractor Ch4 is possibly through a boundary crisis (jHilbornl 
an unstable periodic solution hits the basin boundaries of Ch4- This unstable periodic 
solution possibly has connections with the unstable periodic solutions originating from 
the branch point CB2 of Fig. [5] 

In the next subsection, we will discuss the bifurcations associated with the solutions 
arising from the branch point CB2 whose flow patterns resemble the wavy rolls. 
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Figure 12. Bifurcation diagram Bif-B with the same color convention as BH-13M (Fig. [6}. A 
large chaotic attractor represented by gray dots Ch4 is born at CB3. The chaotic attractor Ch4 
coexists with Chi, Ch2, Ch3 (shown in the inset using brown dots), and SQOR. 




Figure 13. Phase space projection of coexisting SQOR (green) and chaotic attractor Ch4 (gray) 
on the Wioi-Woii plane at r = 1.0342: (a) obtained from model and (b) obtained from DNS. 

4.5. Bifurcation diagram Bif-C of the 27-mode model 

The 27-mode model has another chaotic attractor near the onset that originates from the 
bifurcation point CB2 of Fig. [5] Recall from i M.ll that CB2 is a codimension-2 bifurcation 
point whose stability matrix has a simple zero eigenvalue and an imaginary pair (0, ±iu>). 
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Figure 14. Bifurcation diagram Bif-C of the wavy roll solutions. The blue, red, and brown 
points represent periodic, quasiperiodic and chaotic (Ch5 and ChS) solutions respectively. 

As a consequence, an unstable limit cycle is generated as r is increased beyond 1. The 
attractors from this branch yield periodic, quasiperiodic, and chaotic wavy rolls. See 
Fig- HI for an illustration of a periodic wavy roll. The diagram Bif-C (Fig. [14]) illustrates 
the bifurcation scenario of this type of solutions. The limit cycles generated through 
this bifurcation are unstabl e and they have an unsta b le tor us associated with them as 



discussed in Section 7.4 of 



Guckenheimer fc Holmes! ([1983). As a result, four chaotic 



attractors named ChS are born. A phase space projection of one of the Ch5 attractor is 
shown in Fig. 115( a). Its chaotic nature is ascertained by the broad-band power spectrum 
exhibited in Fig. [TST b). As r is increased further, the size of Ch5 increases till r ~ 1.009 
after which a s ingle large ch aotic attractor Ch6 is generated through an "attractor- 



merging crisis" (jHilborn 



1994D . A phase space projection and power spectrum of Ch6 are 
shown in Fig. fT5T c.d) respectively. The chaotic attractors Ch5 and Ch6 are illustrated in 
the bifurcation diagram Bif-C. 



As r is increased further, we observe a series of phase-locked and stable quasiperiodic 
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Figure 15. Phase space projections on the Wioi-VFiii plane for the wavy rolls and their corre- 
sponding power spectra obtained from the model: (a,b) the chaotic attractor Ch5 at r = 1.005; 
(c,d) the chaotic attractor Ch6 at r = 1.05; (e,f) frequency locked state with /1//2 « 5 at 
r = 1.078; (g,h) quasiperiodic state with /1//2 ~ 6.33 at r — 1.10; (i,j) periodic state at 
r = 1.15 (/l = 4.3). 

solutions. Phase space projections and power spectra of the phase-locked (r = 1.078) and 
quasiperiodic (r = 1.10) states are shown in Fig. ITST e.f) and Fig. fT5Tg.h) respectively. 
The phase-locked and quasiperiodic states are also evident in the bifurcation diagram 
Bif-C as banded and filled regions. Further increase in r transforms the quasiperiodic 
states to a limit cycle through an inverse Niemark-S acker bifurcation. A phase space 
projection and power spectrum of a limit cycle generated for r = 1.15 are shown in 
Fig. [TBTi.i). Note that the transformations of phase-locked, quasiperiodic, and periodic 
states to one another can be understood by the movement of the Floquet multipliers of 
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the underlying limit cycles (see §4.3j) . The above states have also been observed in DNS. 
For example, Fig. [TO] illustrates chaotic, quasiperiodic, and periodic states obtained in 
DNS for r = 1.05, 1.09, and 1.10 respectively. 

When we examine the active modes of Bif-C, we find that only the green colored 
modes of Fig. [T] are active for these convective patterns. The bifurcation diagram Bif-C 
has been generated by setting these modes as nonzero and all other modes as zero. The 
most important among the active modes of Bif-C are (1, 0, 1), (0, 1, 0), and (1, 1, 1) which 
are instrumental for the generation of wavy rolls along the y axis. Naturally, the wavy 
rolls along the x axis will have the complimentary set of modes, e.g., (1,0,0) instead of 
(0,1,0), etc. 

Note that the above sets of solutions (Bif-A, Bif-B, Bif-C) are observed for different sets 
of initial conditions. For both DNS and the model, a random initial condition generally 
produces Bif-B whose basin of attraction appears to be rather large. If we keep only the 
green colored modes of Fig. Q] as nonzero, we obtain Bif-C that correspond to the wavy 
rolls. In fact, Bif-C is generated by first constructing a limit cycle for r — 1.25, and then 
continuing the solution for lower r values. The bifurcation diagram Bif-A is generated 
by starting from the SQ pattern for r = 1.25 for which all the modes except Win, #110, 
#112, #100) #oio, #in, #210, #120, #201, #021, #102, and Z i2 are nonzero. Continuation 
of the above solution however generates ASQ solutions followed by HI, NS1, H2, and 
NS2 bifurcations where all the 27 modes are present. 

Wavy rolls are one of most studied convective patterns in experiments and numerical 
simulations. The bifurcation diagram Bif-C provides a clear explanation for the origin 
of this pattern. In the next section we will provide a quantitative comparison of the 
bifurcation results related to the wavy rolls with those from the experiments and previous 
simulations. 
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Figure 16. Phase space projections on the W101-W111 plane corresponding to the wavy rolls 
obtained from DNS: (a) the chaotic attractor Ch6 at r = 1.05; (b) quasiperiodic state at r = 1.09; 
(c) periodic state at r = 1.10. 

5. Wavy rolls: a quantitative study 

In this section we will analyze time scales of the wavy rolls of Bif-C quantitatively, 
and compare these values with some of the experimental and numerical results. At the 
bifurcation point CB2, the stability matrix has a pair of complex eigenvalues (0, ±iui) 
with uj « 14.2. As a result, the unstable limit cycle originating from CB2 has a time 
period around 2-k/oj w 0.44 in units of d 2 jv (viscous time scale). Subsequent periodic, 
quasiperiodic, and chaotic time series have time scales comparable to the above value 
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since their origin is closely connected to the bifurcation CB2 at r = 1 in Bif-C. Our 
preliminary calculations indicate that the time period of oscillations for these patterns 
are withi n a factor of 10 from this value for r = 1-1.25. 

Earlier iKrishnamurtil (|1973f ) observed time-dependent wavy rolls in her convection ex- 



periments on mercury (P ss 0.02). She observed multiple pe aks with the time p eriod 



ranging from 0.1 to 1 in the time units of d 2 jv (see Fig. 3 of 



Krishnamurti 



1973). Kr- 



ishnamurti's experimental value for the time-scale of th e chaotic wavy rol l s is in the 



same range as our theoretica l time-scale estimated above 



Willis fc Deardorffl dl970h and 



Croquette fc Williams! (j.1989) reported the time period of the oscillatory rolls using their 



experimenta l data f or air (P = 0.7) to be around 1 in the units of d 2 jv. Using numerical 



simulations, 



Lippsl (jl97a) observed ti me periods of osci ll atory rolls to be around 0.24 



0.27 in the units of d 2 / k, for P = 0.7. 



Meneguzzi et al. 



(|1987l ) found the period of the 



wavy oscillations to be around 0.065 viscous time units for P = 0.025. These results are 
in general agreement (within a factor of 10 ) with our theoretical finding based on the 
bifurcation analysis. Note that iBussd (|1972h reported that the time period of the oscil- 



latory instability in the units of d 2 jv is independent of the Prandtl number. Hence the 
time-scales are not expected to vary appreciably even when we change the Prandtl num- 
ber by an order of magnitude which is consistent with the results for mercury (P « 0.02) 
and air (P sa 0.7). Therefore, a comparison of our results for P = with those obtained 
for finite P is also justified. 

Oscillatory instabilities and their saturation through critical bursting have been stud- 



ied by Kumar and coworkers (jKumar. Fauve &: Thual 



1996; 



Kumar. Pal fc Fauve 



2006) 



using several low-dimensional models. They show that the growth of the mode W\q\ is 
saturated by the vorticity mode ^oio- In Fig. [T7Ta.b) we plot the time series of (v\ + v 2 ) 
(sum of kinetic energy along x and z axes) and (v 2 ) (kinetic energy along y axis) com- 
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Figure 17. Time series of (vf +w|)/2 (panel (a)) and {v%) /2 (panel (b)) obtained from the 
model at r = 1.05. Panel (c) shows the time series of the modes Wioi, Win, and Zoio during 
critical bursting. 

puted from our 27-mode model for r = 1.05. These results are in general agreement with 
the results of Kumar and coworkers. The panel (c) of Fig. [17] shows the time series of the 
modes Wioi, Win, and Zqiq that illustrates their growth and subsequent breakdowns 
(critical bursting). Note that the time-scales of oscillations for the modes Wm and Z 010 
are around 0.1, which is in the same range as the theoretical time-scale derived above 
using the bifurcation analysis. 

The above arguments strongly suggest that the origin of the wavy rolls or the oscillatory 
instabilities are intimately related to the purely imaginary pair of eigenvalues at CB2 and 
the limit cycles that originate from it. 
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6. Conclusion 

In conclusion, we explored various flow patterns of zero-P convection and performed a 
detailed bifurcation analysis near the onset using direct numerical simulation and a 27- 
mode low-dimensional model. The low-dimensional model was constructed using the most 
energetic modes of DNS. The results of the DNS and the low-dimensional model are in 
good agreement with each other. Several new chaotic attractors and windows of periodic 
and quasiperiodic rolls have been reported for the first time for zero-P convection. The 
origin and dynamics of all the observed patterns have been explained successfully using 
the bifurcation diagrams. 

The RBC system for P = is chaotic at the onset itself. The stability analysis of 
the 27-mode model indicates three codimension-2 bifurcation points that play critical 
roles in the dynamics of co nvection near th e onset. The chaotic attractors Chi, Ch2, and 
Ch3, described earlier by 



Pal et al. 



(2009), and Ch4 are all related to the bifurcation 
point CB3. Beyond Chi and Ch4, we observe SQOR, OASQ, ASQ, SQ etc., that are 
common to Pal et al.'s 13-mode model. The other codimension-2 bifurcation point, CB2, 
generates chaotic attractors Ch5 and Ch6, and the subsequent periodic, quasiperiodic, 
and phase-locked convective states, which correspond to the wavy roll patterns of RBC 
observed earlier in experiments and simulations. In addition, we find that the frequency 
of the wavy rolls are connected to the imaginary eigenvalues of the stability matrix at 
the CB2 bifurcation point. Thus, the bifurcation analysis presented in the paper provides 
useful insights into the origin of the wavy rolls of RBC. 



Inte restingly, the bifurcation diagram of 30-mode model of 



Mishra. Wahi fc Verma 



( 20101 ) for P = 0.0002 matches quite closely with Bif-A of our model. This reinforces 



earlier observations that zero-P convection is a valid limit of low-P convection as P — » 



(Thual 



1992). The extension of the present study to low-P convection in relation to wavy 
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rolls will be very valuable for understanding the experimental and numerical findings 

near the onset. 

We are thankful to Krishna Kumar and Pankaj K. Mishra for useful discussions. This 
work is supported by the Swarnajayanti fellowship grant to MKV by Department of 
Science and Technology, India. 
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